% lapse+edit: edited to allow for endo lapses

function [v, lapse_NN, LL] = predict_v_lapse(p1, p2_input, ...
    alpha,  Delta_xi,...
    N_y, resource_y1, resource_y2, pr_s_row, KM,...
    rho, crra, beta_c, U0_mm, B1_c, B2_c, T1, delta)

NN = size(p1,1); % major insurers + # fringes

v = zeros(N_y, NN); 
c1_all = zeros(N_y, NN);

lapse_NN = zeros(N_y, NN);% lapse probability with k intergrated over

% loop thru each insurer jj in the market
for jj=1:NN
    
    % 2nd period probability distribution in state_mm -> expanded 
    pr_s = pr_s_row'; % KM x 1

    % jj's prices
    p_jj1 = p1(jj);
    p_jj2 = p2_input(jj, :); 
    
    Delta_xi_jj = Delta_xi(jj);
    
    % first-period annual consumption (N_y x 1)
    c1 = resource_y1-p_jj1;

    % second-period annual consumption in each state k if not lapse (N_y x KM)
    c2 = repmat(resource_y2,1,KM)-repmat(p_jj2,N_y,1);

    % second-period utility in each state k if not lapse (N_y x KM)
    U_stay = B2_c*(alpha*u_c(c2,rho,crra) +  Delta_xi_jj);

    % cannot afford p2
    force_lapse = (c2<1);

    % second-period utility in each state k if lapse 
    U_lapse = U0_mm; % N_y x K
    U_lapse_mat = repmat(U_lapse, 1, 4); % N_y x KM

    % expected utility conditional on realized state 
    U = log(exp(U_stay) + exp(U_lapse_mat));% N_y x KM
    U(force_lapse==1) = U_lapse_mat(force_lapse==1);

    % lapse+edit: predict lapse probability
    lapse_mat = exp(U_lapse_mat)./(exp(U_stay) + exp(U_lapse_mat)); %N_y x KM
    lapse_mat(force_lapse==1) = 1;

    % save N_y x KM lapse prob for insurer jj in structure
    LL(jj).lapse_mat = lapse_mat; % N_y x KM

    % - integrate over 2nd-period states
    lapse_vec = lapse_mat*pr_s; % N_y x 1
    lapse_NN(:,jj)=lapse_vec;

    % expected 2nd-period utility (N_y x 1)
    EU = U*pr_s;
    
    % size(v_jj) = N_y x 1
    v_jj = B1_c*(alpha*u_c(c1,rho,crra) +  Delta_xi_jj) + (beta_c^T1)*EU;

    % set v_jj to -Inf if income less than p_jj in any state
    flag = [c1 c2];
    flag = (flag<0);
    flag = (sum(flag,2)>0);
    v_jj(flag==1)=-Inf;
    
    % save
    v(:,jj) = v_jj;
    c1_all(:,jj) = c1;
end


